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Introduction 


The  objectives  of  this  3 -year  project  were  to  formulate  high  fidelity  computational  models  for 
soft  biological  tissues  that  can  eventually  be  implemented  within  existing  software  for  use  in 
telemedicine  and  surgical  simulation.  We  proposed  to  implement  quasi-linear  viscoelastic  theory 
(QLV)  in  3D,  develop  fractional-order  viscoelasticity  (FOV)  to  a  higher  level,  and  implement 
more  sophisticated  3-D  models  into  simulations  of  soft  tissue  behavior.  As  a  laboratory  with  long 
experience  in  experimental  work  and  validation,  we  also  proposed  to  perform  bench  studies  and 
tissue  loading  experiments  so  that  these  models  could  be  properly  validated. 


Body 

Our  general  approach  can  categorized  into:  (i)  theoretical  development,  (ii)  algorithmic 
development/coding,  (iii)  code  validation  through  analytical  solutions  when  these  exist,  (iv) 
theory/code  validation  through  comparison  with  experimental  data  and,  (v)  simulating  a  number 
of  simple  but  realistic  surgical  processes.  In  this  report,  the  individual  accomplishments  will  be 
reported  in  separate  sections. 


(a)  Improved  parameter  estimation  methods  and  novel  experimental  techniques 


A  new  parameter  estimation  method  for  QLV  that  is 
less  sensitive  to  imperfections  in  the  experimental 
data  has  been  developed  and  published  '.  The 
method  is  robust,  insensitive  to  the  loading  history 
of  the  test  specimen  and  handles  noise  in  the 
position  channel  quite  well.  Although  this  is  a  one¬ 
dimensional  approach,  it  is  remarkably  predictive  of 
the  time-varying,  one  dimensional  response  of 
elongated  structures,  such  as  ligaments,  tendons,  and 
strips  of  tissue.  This  approach  would  thus  be  useful 
for  modeling  surgery  of  tissue  strips  with  the 
appropriate  geometry. 

(b)  New  theoretical  and  numerical  developments  for 
fractional-order  viscoelasticity  (FOV) 

Most  previous  studies  of  soft  tissue  viscoelasticity 
have  used  Boltzmann’s  general  theory  of 
viscoelasticity,  extended  to  include  nonlinear  elastic 
behavior  and  a  continuous  box-spectrum  relaxation 
response  (Fung’s  well-known  quasilinear 
viscoelasticity  theory,  QLV).  We  have  introduce  an 


Time  (seconds) 


Time  (seconds) 


Figure  1:  Plots  of  the  cyclic  response  data  and  model 
predictions  for  a  typical  specimen  for  (a)  all  20  cycles, 
and  expanded  to  the  first  3  cycles  (b&c).  Also  shown 
is  the  pointwise  RMS  error  for  both  methods  (d). 


1  Doehring  T.C.,  Carew  E.O.,  Vesely  I.  The  effect  of  strain  rate  on  the  viscoelastic  response  of  aortic  valve  tissue:  a 
direct-fit  approach.  Ann.Biomed.Eng.  32(2):223-32,  2004. 
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alternative:  the  fractional  order  viscoelasticity  (FOV)  theory,  which  uses  a  fractional  integral  for 
the  relaxation  response.  FOV  has  some  appealing  characteristics,  such  as  more  flexible  control 
over  the  early  relaxation  behavior.  FOV  also  implies  a  hierarchical  (fractal-like)  structure,  which 
is  relevant  to  connective  tissues2 *.  The  theoretical  work  regarding  the  implementation  of  FOV  has 
been  completed  by  our  German  colleague,  Dr.  Kai  Diethehn  and  two  publications  have  resulted 
3,4.  We  validated  the  fidelity  of  this  method  by  comparing  a  predicted  viscoelastic  response  of  a 
modeled  experiment  with  the  real  acquired  data  and  found  errors  typically  below  2%  (see  Fig.  1). 
Overall,  FOV  is  only  slightly  better  than  QLV  in  predicting  long  term  stress  relaxation  and  cyclic 
loading,  but  shows  promise  because  of  its  hierarchical  nature  that  is  more  representative  of  the 
internal  complexity  of  biological  tissues. 

(  c)  Validation  of  the  Anisotropic  Hyperviscoelastic  (ANHVE)  model  with  soft  biological  tissues 

Early  in  the  program,  we  contracted  Dr.  Saleeb  at  the  University  of  Akron  to  implement  his 
ANHVE  model  (previously  referred  to  as  the  LSVEA  formulation)  to  soft  biological  tissues. 
This  model,  by  definition,  is  a  phenomenological  model,  in  that  the  various  parameters  estimated 
for  the  multi-mechanism  model  of  the  material  in 
question,  do  not  have  any  physical  representation. 

We  nevertheless  pushed  ahead  with  this  approach 
because  it  is  computationally  expedient. 

Validation  of  the  ANHVE  model  involved 
parameter  estimation  from  biaxial  (2D) 
experiments  and  subsequent  simulations  of 
independent  uniaxial,  biaxial  (planar)  and  inflation 
(nonplanar)  tests.  The  response  predicted  by  the 
numerical  approach  was  reasonable  (see  Fig.  2)  but 
refinement  is  needed.  Tuning  of  the  model  for 
more  complex  structures  became  difficult,  so  we 
thus  opted  to  develop  our  own  model, 
incorporating  features  of  the  tissue  that  are  more 
directly  observable  -  the  fiber  distribution 
population  -  as  shown  in  the  next  section. 

(d)  Invariant  theory  for  dispersed  transverse  isotropy. 

The  past  few  years  have  seen  an  increased  interest  in  nonlinear  continuum  mechanics  as  a 
framework  for  describing  the  mechanical  behavior  of  soft  biological  tissues.  The  mathematics  in 
this  field  are  now  well-established  and  thus  provide  a  foundation  for  thennodynamically- 
consistent  constitutive  equations.  These  equations  are  readily  handled  by  the  theory  of 
invariants,  and  incorporate  finite  deformations,  geometric  and  material  nonlinearities  and  tissue 
anisotropy.  Powerful  computers  and  well-developed  computational  techniques  now  make 

2  Doehring  T.C.,  Freed  A.H.,  Carew  E.O.,  Vesely  I.,  Fractional  order  viscoelasticity  of  the  aortic  valve:  An 
alternative  to  QLV.  J.  Biomech.  Eng.  127(4):700-8,  2005. 

’  Diethelm  K.,  Ford  N.J.,  Freed  A.D.,  Luchko  Yu.,  Algorithms  for  the  fractional  calculus:  A  selection  of  numerical 
methods.  Computer  Methods  in  Applied  Mechanics  and  Engineering,  194:734-773,  2005. 

4  Diethelm  K.,  Weilbeer  M.,  A  numerical  approach  for  Joulin’s  model  of  appoint  source  initiated  flame.  Fractional 
Calculus  and  Applied  Analysis,  (7):2,  191-212,  2004. 


Figure  2:  Image  (left)  of  a  presurized  segment  of  aortic 
valve  tissue  showing  how  the  markers  on  the  surface 
more  apart  in  response  to  the  distension  of  the 
membrane.  An  FEA  representation  of  the  same 
experiment  (right)  showing  the  anisotripic  deformation 
of  the  tissue  segment,  being  different  in  the  two 
principle  fiber  directions. 
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simulations  of  tissue-level  mechanics  feasible.  Interactions  between  crossing  fiber  populations 
have  been  implicated  in  generating  the  complex  biaxial  material  properties  of  most  biological 
tissues.  For  example,  under  certain  combinations  of  biaxial  tension,  the  material  can  shrink  in 
one  principle  direction  of  applied  stress  Our  recently  developed  a  constitutive  equation  for  soft- 
tissues  based  on  a  new  invariant  theory  for  dispersed  transverse  isotropy5.  The  current  model 
replaces  a  numerically  evaluated  integral  with  a  single  analytic  scalar,  and  is  thus 
computationally  expedient,  increasing  the  speed  of  our  simulations  three  fold. 

In  our  first  validation,  we 
made  use  of  our  long  standing 
expertise  in  heart  valve 
biomechancis,  and  simulated 
the  opening  and  closing 
behavior  of  an  artificial  heart 
valve.  The  first  part  of  the 
simulation  involved  first 
fitting  the  model  to 
experimental  data  of  a  series 
of  biaxial  tensile  tests  of  the 
leaflet  material.  The  fit  of  our 
model  to  these  experimental  data  (see  Fig.  3)  demonstrates  that  our  model  can  simulate  both  the 
non-linear  stress/strain  behavior  of  this  tissue,  as  well  as  the  lateral  contraction  that  occurs  under 
certain  combinations  of  biaxial  loading.  The  resulted  simulation  of  making  use  of  this  material 
model  showed  very  good  agreement  to  know  behavior  of  heart  valve  tissues. 

(e)  Micromechanical  model 

To  obtain  the  greatest  fidelity  of  a  model  to  the  real  material,  one  must  take  into  account  the 
internal  microstructure.  To  that  end,  we  have  been  working  with  Drs.  Pindera  and  Aboudi,  from 
the  University  of  Virginia  to  complete  the  development  of  a  high-fidelity  micromechanical 
model.  This  model  uses  the  principle  of  Generalized  Method  of  Cells  (GMC)  and  can  predict 
both  the  macro-level  response  and  the  micro-level  stress  and  deformation  fields  in  the  individual 
material  phases.  GMC  assumes  that  a  region  of  material  space  consists  of  a  periodic  array  of 
cells  with  identical  microstructure  and  material  properties.  A  "homogenization  technique"  is  used 
to  generate  equations  for  the  mechanical  properties  of  each  cell  that  result  from  its  complex 
microstructure.  These  equations  represent  a  new  material  for  which  the  aggregate  properties 
evolve  as  the  microstructure  changes  with  the  application  of  strain.  The  micromechanical  model 
can  run  in  parallel  with  the  main  Finite  Element  Analysis  (FEA)  code,  generating  what  is 
essentially  a  new  stiffness  matrix  for  the  material  at  each  time  step.  This  approach  enables  a 
complex  fibrous  material  to  be  represented  by  relatively  few  elements,  thus  improving 
computational  performance  without  loss  of  mechanical  fidelity. 

To  validate  this  model  and  to  explore  its  utility  in  simulating  the  complex  micromechanical 
behavior  of  biological  tissues,  we  have  again  turned  to  heart  valve  tissue  as  the  object  of  interest. 


Figure  3:  Simulations  of  standard  biaxial  tissue  tests  produce  a  family  of 
stress/strain  curves  that  correspond  closely  to  the  raw  data  (left).  Finite 
element  simulations  based  on  the  new  model  demonstrate  physiologic 
patterns  of  stress  and  leaflet  configuration  in  open  and  closed  position. 


5  Freed  A.D.,  Einstein  D.R.,  Vesely  I.,  Invariant  formulation  for  dispersed  transverse  isotropy  in  aortic  heart  valves: 
An  efficient  means  for  modeling  fiber  splay.  Biomechanics  and  Modeling  in  Mechanobiology,  4:100-1 17,  2005. 
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We  have  shown  before  that  the  modulus  of  Marginal 

mitral  valve  chordae  varies  with  chordal  A 
type.  Even  though  the  total  amount  of 
collagen  is  the  same  for  all  three  types,  the 
fibril  size  distribution  is  different  (Fig. 

4A,B).  In  a  preliminary  study  to  test  the 
predictive  ability  of  a  micromechanical 
analysis,  we  simulated  both  the  effects  of 
collagen  crimp  and  collagen  fibril  diameter. 

Collagen  crimp  is  expected  to  affect  the 
ultimate  extensibility  of  these  tissues  and 
collagen  fibril  distribution  is  expected  to 
affect  the  tangent  modulus,  as  predicted  by 
our  structural  model.  Although  these  two 
phenomena  occur  at  different  scales,  for  the 
sake  of  simplicity,  we  generated  structural 
models  that  incorporated  both  at  once.  The 
collagen  crimp  was  set  to  16.8,  12.6  and 
10.8  mm  for  basal,  marginal  and  strut 
chordae  (Fig.  4C).  The  collagen  fibril  size 
distribution  was  composed  of  various  ratios 
of  thick  and  thin  fibers  of  diameter  2,  4  and 
6.  Different  combinations  of  these  fibers 
were  placed  in  the  simulation  grid  (Fig. 

4C),  to  mimic  the  fibril  size  distribution 
shown  in  Fig.  4B,  albeit  at  a  much  coarser 
resolution.  In  each  case,  the  volume 
fraction  of  fiber  and  matrix  was  kept  at 
50%,  so  that  any  mechanical  differences 
resulted  only  from  the  differing 
microstructure.  Starting  values  for 
collagen  fibril  stiffness  were  inferred  from 
force  spectroscopy  data  while  matrix 
properties  were  estimated  from  pre¬ 
transition  stress-strain  data.  The  model  was 
loaded  incrementally  and  an  effective 
stress/strain  curve  for  the  “material”  was 
generated.  Interestingly,  this  very  simple 

model  demonstrated  the  same  pattern  of  increasing  extensibility  and  decreasing  tangent  modulus 
in  the  order  marginal,  basal,  strut  (Fig.  4D),  exactly  as  was  measured  experimentally.  Probing 
deeper  into  the  stress  environment  between  the  fibrils  (Fig.  4E)  reveals  the  reason  why  chordae 
with  a  lower  crimp  period  (more  crimped  collagen)  and  larger  diameter  collagen  fibrils  have  a 
lower  tangent  modulus.  The  larger  fibrils  bear  a  greater  proportion  of  the  stress  than  the  thinner 
fibrils  at  a  lower  extensions  -  essentially  when  they  are  still  more  crimped  and  thus  softer.  This 
kind  of  insight  would  have  been  difficult  to  obtain  without  access  to  micromechanical  models. 


10.0 
Strain  (%) 

Figure  4:  (A)  Electron  microscopy  images  of  collagen  fibrils 
in  marginal,  basal  and  strut  chordae  and  their  corresponding 
fibril  diameter  distributions  (B).  (C)  The  corresponding 
micromechanics  model  of  each  chordae.  (D)The  model 
demonstrated  the  same  reduction  in  stiffness  in  the  order  of 
marginal,  basal,  and  strut,  as  did  the  ohginal  mechanical 
data  (E)  The  distribution  of  the  internal  stresses,  at  a  total 
stress  of  2MPa,  shows  that  the  thin  fibrils  (arrows)  in  the 
strut  chordae  are  stressed  less  than  the  thin  fibrils  in  the 
marginal  chordae,  pointing  to  possible  mechanism  for  the 
stiffness  variation. 
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(f)  Implementation  of  inverse  methods 

The  advantage  of  a  constitutive  model  is  speed,  at  the  expense  of  fidelity.  The  advantage  of  a 
micromechanics  approach  is  fidelity,  possibly  at  the  expense  of  speed.  To  obtain  the  best  of  both 
approaches,  we  have  considered  hybrid  approaches.  A  hybrid  approach  between  constitutive  and 
microstructural  modeling  is  regional  tuning  of  material  properties  using  inverse  methods. 
Connective  tissues,  such  as  heart  valves  and  blood  vessels,  have  spatially  varying  material 
properties.  However,  a  small  set  of  constitutive  models  can  be  used  for  the  entire  structure  if 
their  constants  are  tuned  so  that  the  mechanics  of  the  individual  finite  elements  represent  the 
mechanics  of  tissue  at  the  same  location.  In  preparation  for  this  project,  we  have  developed  an 
inverse  method  based  on  the  response  surface  methodology  (RSM)  that  can  be  coupled  to  whole 
tissue  loading  experiments6. 

Like  all  inverse  FEA  methods,  RSM  is  driven  by  a  series  of  forward  FEA  simulations  that 
generate  synthetic  data  whose  deviation  from  the  real  experimental  data  needs  to  be  minimized. 
The  solution  space  has  as  many  dimensions  as  there  are  parameters  to  be  tuned,  and  the  response 
surface  usually  contains  many  local  minima.  The  response  surface  methodology  first  constructs 
a  low  order  approximation  to  the  response  from  a  small  number  of  forward  simulations 
distributed  through  the  solution  space.  The  minimum  of  one  response  surface  becomes  the 
starting  point  for  the  next  set  of  simulations.  The  parameter  space  is  contracted  with  each 
iteration  so  that,  in  the  neighborhood  of  the  global  minimum,  finer  features  of  the  response  are 
captured.  This  cycle  of  forward  solutions  followed 
by  the  minimization  of  the  objective  function  is 
repeated  until  some  minimal  residual  is  reached. 

The  approximate  nature  of  the  response  surface  over 
a  wider  parameter  space  at  the  beginning  makes  it 
relatively  insensitive  to  local  minima  and 
experimental  noise.  Because  minimization  proceeds 
along  the  approximate  response  surface,  RSM  is  an 
inherently  parallel  method,  and  thus  is  well-suited  to 
massively  parallel  computer  architectures.  The 
actual  tuning  can  be  done  on  data  sets  from  isolated 
tissues,  whole  vessel  loading  experiments  or  3-D 
image  data  from  patients.  Indeed,  in  the  experiment 
shown  in  Fig.  5,  we  see  that  a  simulation  of  vessel 
clamping  (one  of  the  objectives  of  this  project), 
shows  a  pattern  of  curvature  in  the  clamped  zone 
that  is  very  similar  to  the  curvature  of  the  real  tissue 
being  clamped  (in  this  case  a  pig  aorta).  Where 
inverse  methods  are  particularly  useful  is  the 
blending  of  experimental  data  from  many  sources 
into  a  single  material  model.  For  example,  the 
pressurization  experiments  shown  in  Fig.  2  can  be 
used  as  target  data  along  with  the  experimental  data 


6  Einstein  D.R.,  Freed  A.D.,  Stander  N.,  Fata  B.,  Vesely  I.,  Inverse  Parameter  Fitting  of  Biological  Tissues:  A 
response  surface  approach.  Ann.Biomed.Eng.,  33(12):1819-1830,  2005. 


Figure  5:  Image  (bottom)  of  a  presurized  segment  of 
pig  aorta  being  clamped  with  forceps,  simulating  one 
of  the  key  features  of  virtual  surgery  training  - 
manipulating  tissues.  The  image  at  the  top  is  an  FEA 
simulation  of  the  same  condition,  complete  with 
fluid/solid  interactions.  By  matching  the  deformed 
profiles  of  the  FEA  simulation  to  real  experimental 
data  of  the  same  deformation,  confidence  can  be  had 
that  the  model  matches  reality.  Tuning  of  such  models 
can  be  done  best  through  the  inverse  FEA  methods. 


of  Fig.  5,  to  arrive  at  a  material  model  that  can  satisfy  both  conditions. 

Key  Research  Accomplishments 

Computational  modeling  of  soft  biological  tissues  is  extremely  challenging.  Unlike  conventional 
materials,  where  the  stress/strain  response  is  well  characterized,  biological  materials  are  poorly 
understood.  In  the  absence  of  a  good  experimental  knowledge  of  tissue  mechanics,  formulation 
of  constitutive  models  becomes  even  more  challenging.  Coupled  with  this  lack  of  knowledge 
over  detailed  mechanics  is  a  poor  understanding  of  microstructure.  Without  an  understanding  of 
microstructure,  complexity  cannot  be  build  into  the  models  to  improve  predictive  fidelity. 
Progress  in  this  field  has  therefore  come  about  by  an  iterative  improvement  in  both  modeling 
capability  and  an  understanding  of  tissue  structure  and  mechanics.  We  have  become  the  leaders 
in  this  field,  largely  because  of  the  support  form  the  US  Army.  Over  the  3  years  of  support 
(which  has  been  extended  due  to  the  disruption  of  activities  by  the  move  from  Ohio  to 
California),  we  have  developed  the  following  key  pieces  of  technology: 

(i)  A  one-dimensional  model  of  fractional  order  viscoelasticity  that  is  representative  of 
the  hierarchical  nature  of  complex  biological  tissues.  This  is  the  first-of-its-kind 
implementation  of  fractional  order  calculus  in  the  analysis  of  biomaterials. 

(ii)  A  micromechanical  approach  to  modeling  soft  biological  tissues  that  incorporates 
material  and  geometrical  non-linearity.  This  is  also  a  first-of-its-kind  application  of 
GMC  to  biological  materials. 

(iii)  A  computationally  expedient  and  very  accurate  constitutive  representation  of  the 
dispersed  isotropy  that  is  typical  of  biological  tissues.  This  method  is  faster  than  other 
methods  currently  being  used. 

These  three  key  pieces  of  technology  form  the  foundation  with  which  to  study  soft  tissue 
behavior  in  more  detail  into  the  future. 

Reportable  Outcomes 

The  reportable  outcomes  of  this  research  are  the  publications  listed  in  the  Footnotes  of  this 
report.  These  publications  are  summarized  in  the  References  section  below.  We  believe  that  we 
have  satisfied  the  main  objectives  of  this  research  program  -  to  generate  a  collection  of  novel, 
more  accurate  computational  models  and  to  validate  them  using  a  series  of  materials  tests, 
including  a  vascular  clamping  experiement. 

Conclusions 

These  models,  however,  are  not  yet  ready  for  implementation  into  surgical  simulators.  They  are 
still  much  too  slow  for  real  time  use.  We  expect  that  in  the  future,  we  and  others  will  take  these 
models  and  implement  them  in  faster  hardware  and  software,  and  also  simplify  them  so  that  the 
essence  of  the  features  of  soft  tissues  that  are  necessary  for  their  use  in  surgical  simulation,  can 
be  finally  implemented  in  real  training  systems. 
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